N2O AF fitting tests (data analysis plots only)¶

11/04/24

WAS WORKING FOR BASIC EXPORT AFTER COPY, BUT NOW FAILING - probably due to running kernel... NEED TO ADD SETUP STUFF BACK IN...

Notes¶

Results only from test fitting, 1000 fits + noise using 10K ADMs, 41 t-point over half-revival, and adding ~10% noise.

General comments:

  • BLM(t) fitting very good, very small spread in best fit outputs.
  • Some magnitudes and phases retrieved, but generally poor agreement with inputs.
  • Easiest way to visulalise is via density matrices: shows good results for some values, particularly for imaginary terms.
    • Poor results for many real (magnitude) terms, esp. odd-L, likely because this is not defined in symmetrized dataset.
    • The main feature in the theory results around L=2,0 looks OK, but then get similar mag terms for (1,0), (3,0) and (4,0) terms, appearing as grid pattern in recon results.
    • TODO: try symmetrizing theory results to compare properly?
    • TODO: try allowing large L terms to fall off (e.g. Wigner threshold law) in fitting.

  • Data set: /home/jovyan/N2O/epsman2024/proc-fock/N2O-fitTest_1000_080424/N2O_559_fit_orb9_8.1_noise_080424_09-45-35.pickle
  • Running code (with fits): N2O/epsman2024/proc-fock/N2O_AF_channelFuncs_240324-Fock-fitting.ipynb
In [1]:
# Quick hack to override default HTML template
# NOT required in some JLab versions.
# https://stackoverflow.com/a/63777508
# https://stackoverflow.com/questions/21971449/how-do-i-increase-the-cell-width-of-the-jupyter-ipython-notebook-in-my-browser
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
/tmp/ipykernel_48708/28059362.py:5: DeprecationWarning: Importing display from IPython.core.display is deprecated since IPython 7.14, please import from IPython display
  from IPython.core.display import display, HTML
In [2]:
# %matplotlib inline
# import matplotlib.pyplot as plt
In [3]:
# %matplotlib notebook
# Throws IPython errors?
In [4]:
# import matplotlib.pyplot as plt
# plt.style.use('classic')
# import numpy as np

# %matplotlib inline

Load alignment data¶

Now with ADM class.

In [5]:
# from epsproc.plot import hvPlotters
# hvPlotters.setPlotDefaults(fSize = [700,300], imgSize = 600)
In [6]:
from epsproc.classes.alignment import ADM
from pathlib import Path

# fileBase = Path('/home/jovyan/jake-home/ePS/N2O/N2O_valence') # Data dir on Jake (from ePSprc container)
fileBase = Path('/home/jovyan/N2O/N2O_valence') # Data dir on Jake (from ePSprc container)
# dataPath = Path('~/ePS/N2O/alignment/Wigner D polyatomic normalized')
ADMdataPath = Path('/home/jovyan/N2O/alignment/Wigner D polyatomic normalized')
ADMtype = '.txt'

ADMin = ADM(fileBase=ADMdataPath.expanduser(), ext = ADMtype)

ADMin.loadData(keyType='setADMsT', normType='wignerDpoly')
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
* sparse not found, sparse matrix forms not available. 
* natsort not found, some sorting functions not available. 
* Setting plotter defaults with epsproc.basicPlotters.setPlotters(). Run directly to modify, or change options in local env.
* Set Holoviews with bokeh.
* pyevtk not found, VTK export not available. 
Scanning /home/jovyan/N2O/alignment/Wigner D polyatomic normalized for ePS jobs.
Found ePS output files in subdirs: [PosixPath('/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A400'), PosixPath('/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A800'), PosixPath('/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A600'), PosixPath('/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A200')]

*** Scanning dir
/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A400
Found 6 .txt file(s)


*** Scanning dir
/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A800
Found 6 .txt file(s)


*** Scanning dir
/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A600
Found 6 .txt file(s)


*** Scanning dir
/home/jovyan/N2O/alignment/Wigner D polyatomic normalized/A200
Found 6 .txt file(s)

Set self.norm from self.norms['wignerDpoly'].
Found t-index /home/jovyan/N2O/alignment/time.txt
In [7]:
ADMin.plot(width=1000)
# ADMin.plot(xlim=(30,50), width=800)  # Auto-stack plus pass args & zoom in on specific region (note slider will reset region with interactive zoom)

Load matrix elements¶

In [8]:
# Plotters
# from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob

# For fitting/basis exploration use PEMtk base class...
from pemtk.fit.fitClass import pemtkFit

import warnings
# warnings.filterwarnings('once')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore')   # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [9]:
# # Scan for subdirs, based on existing routine in getFiles()

# fileBase = Path('/home/paul/ePS/OCS/OCS_survey')  # Data dir on Stimpy

# fileBase = Path('~/fock-mount/globalhome/eps/N2O/N2O_valence').expanduser() # Data dir on Jake
In [10]:
# TODO: fix orb label here, currently relies on (different) fixed format

# ePSproc - read all files
# data = ePSmultiJob(fileBase, verbose = 0)
# data.scanFiles()

# For PEMtk case - currently need to set specific data file
data = pemtkFit(fileBase, verbose = 0)
data.scanFiles(fileIn=[Path(fileBase,'orb8_S',
                            'N2O_valence.orb8_S_E0.1_4.0_28.1eV.out').expanduser().as_posix()])
data.jobsSummary()
# TODO: fix labelling here - have default +1 added?
*** Job subset details
Key: subset
No 'job' info set for self.data[subset].

*** Job orb9 details
Key: orb9
Dir /home/jovyan/N2O/N2O_valence/orb8_S, 1 file(s).
{   'batch': 'ePS N2O, batch N2O_valence, orbital orb8_S',
    'event': 'orb 8 (S/CAv) ionization, batch N2O_valence, None.',
    'orbE': -19.00715329282262,
    'orbLabel': 'S/CAv'}

Compute AF¶

Note default case assumes linear pol, $z$-axis alignment.

In [11]:
# Set ADMs to use for calc - set for specific T and downsample

# Subselection, note temperature is set by dataKey here
temp = '10K'
ADMin.subsetADMs(dataKey = temp, trange=[35,45],tStep = 5)

# Plot subselection
ADMin.plot(keys='ADM')

# Set to main data structure for calcs
data.data['ADM'] = ADMin.data['ADM']
Selecting 41 points
In [12]:
data.data['orb9']['matE'].Eke
Out[12]:
<xarray.DataArray 'Eke' (Eke: 8)>
array([ 0.1,  4.1,  8.1, 12.1, 16.1, 20.1, 24.1, 28.1])
Coordinates:
  * Eke      (Eke) float64 0.1 4.1 8.1 12.1 16.1 20.1 24.1 28.1
    Ehv      (Eke) float64 19.1 23.1 27.1 31.1 35.1 39.1 43.1 47.1
    SF       (Eke) complex128 (2.3800419+3.3895553j) ... (3.7370686+2.1587197j)
Attributes:
    units:    eV
xarray.DataArray
'Eke'
  • Eke: 8
  • 0.1 4.1 8.1 12.1 16.1 20.1 24.1 28.1
    array([ 0.1,  4.1,  8.1, 12.1, 16.1, 20.1, 24.1, 28.1])
    • Eke
      (Eke)
      float64
      0.1 4.1 8.1 12.1 ... 20.1 24.1 28.1
      units :
      eV
      array([ 0.1,  4.1,  8.1, 12.1, 16.1, 20.1, 24.1, 28.1])
    • Ehv
      (Eke)
      float64
      19.1 23.1 27.1 ... 39.1 43.1 47.1
      array([19.1, 23.1, 27.1, 31.1, 35.1, 39.1, 43.1, 47.1])
    • SF
      (Eke)
      complex128
      (2.3800419+3.3895553j) ... (3.73...
      array([2.3800419+3.3895553j, 2.6173399+3.0822452j, 2.8348434+2.8457598j,
             3.0368086+2.6565005j, 3.2261549+2.5005878j, 3.4049881+2.3692545j,
             3.5748864+2.2566545j, 3.7370686+2.1587197j])
  • units :
    eV
In [13]:
# Configuration for data seletion (PEMtk case)

# Matrix element sub-selection
orb = 'orb9'
Eke = 8.1
data.selOpts['matE'] = {'thres': 0.01, 'inds': {'Type':'L', 'Eke':Eke},'sq':False}
data.setSubset(dataKey = orb, dataType = 'matE')  # Subselect from 'orb5' dataset, matrix elements


# # And for the polarisation geometries...
data.setPolGeoms()
data.selOpts['pol'] = {'inds': {'Labels': 'z'}}
data.setSubset(dataKey = 'pol', dataType = 'pol')

# # And for the ADMs...
# data.selOpts['ADM'] = {}   #{'thres': 0.01, 'inds': {'Type':'L', 'Eke':1.1}}
# tRange = [data.data['ADM']['ADM'].t[0], data.data['ADM']['ADM'].t[60], 4]  # Set first 0:60 in steps
# # data.setSubset(dataKey = 'ADM', dataType = 'ADM', sliceParams = {'t':[4, 5, 4]})
# data.setSubset(dataKey = 'ADM', dataType = 'ADM', sliceParams = {'t':tRange})

data.selOpts['ADM'] = {} 
data.setSubset(dataKey = 'ADM', dataType = 'ADM') #, sliceParams = {'t':[220, 580, 5]})
Subselected from dataset 'orb9', dataType 'matE': 114 from 2976 points (3.83%)
Subselected from dataset 'pol', dataType 'pol': 1 from 3 points (33.33%)
Subselected from dataset 'ADM', dataType 'ADM': 205 from 205 points (100.00%)
In [14]:
# PEMtk setup fit + basis set return
import epsproc as ep
from epsproc.geomFunc.geomCalc import setPhaseConventions
phaseCons = setPhaseConventions(phaseConvention = 'E')

BetaNormX, basis = data.afblmMatEfit(phaseConvention=phaseCons)  # Note this returns product basis only.

BetaNormX2, basisFull = ep.geomFunc.afblmXprod(data.data[data.subKey]['matE'], AKQS=data.data[data.subKey]['ADM'],
                                               basisReturn = 'Full', 
                                               thres=None, selDims={}, sqThres=False,
                                               phaseConvention=phaseCons)
In [15]:
%matplotlib inline

import epsproc as ep

ep.lmPlot((basisFull['BLMtableResort'] * basisFull['polProd']), xDim='t', cmap='vlag', thres=1e-3);
Set dataType (No dataType)
Plotting data (No filename), pType=a, thres=0.001, with Seaborn

Channel functions¶

Code from QM3, https://phockett.github.io/Quantum-Metrology-with-Photoelectrons-Vol3/part1/theory_info_content_200723.html#information-content-from-channel-functions

In [16]:
# Define a set of channel functions to test
channelFuncs = (basisFull['BLMtableResort'] * basisFull['polProd'])

# For illustrative purposes, define a subset to use for analysis
channelFuncsSubset = channelFuncs.sel(Labels='A').sel({'S-Rp':0,'mu':0,'mup':0})  #.sel(L=2)

# Check dimensions
print(f"Available dimensions: {channelFuncs.dims}")
print(f"Subset dimensions: {channelFuncsSubset.dims}")
Available dimensions: ('m', 'mp', 'S-Rp', 'l', 'lp', 'L', 'mu', 'mup', 'Labels', 'M', 't')
Subset dimensions: ('m', 'mp', 'l', 'lp', 'L', 'M', 't')
In [17]:
# Convert to PD and tabulate with epsproc functionality
# Note restack along 't' dimension
channelFuncsSubsetPD, _ = ep.util.multiDimXrToPD(channelFuncsSubset.squeeze().real, 
                                                 thres=1e-4, colDims='t')

# Round values to 1 d.p., then apply statistical methods
# Compute per basis index and display
channelFuncsSubsetPD.T.round(2).agg(['min','max','var','count','nunique']).T  
Out[17]:
min max var count nunique
L l lp m mp
0 0 0 0 0 0.0 0.0 0.0 41.0 1.0
1 1 -1 1 0.0 0.0 0.0 41.0 1.0
0 0 0.0 0.0 0.0 41.0 1.0
1 -1 0.0 0.0 0.0 41.0 1.0
2 2 -1 1 0.0 0.0 0.0 41.0 1.0
... ... ... ... ... ... ... ... ... ...
10 5 5 0 0 0.0 0.0 0.0 2.0 1.0
1 -1 -0.0 -0.0 0.0 2.0 1.0
6 4 -1 1 -0.0 -0.0 0.0 2.0 1.0
0 0 0.0 0.0 0.0 2.0 1.0
1 -1 -0.0 -0.0 0.0 2.0 1.0

182 rows × 5 columns

In [18]:
%matplotlib inline  
# NEED TO CALL inline AGAIN HERE, not sure what is going on - something odd in env/Jupyter version?

# channelFuncsSubsetPD.transform(demean,axis='columns') 
# cmap=None   # cmap = None for default. 'vlag' good?
cmap = 'vlag'

# Channel funcs
ep.lmPlot(channelFuncsSubset, xDim='t', cmap='vlag', thres=1e-3);

# Define demean function and apply (from https://stackoverflow.com/a/26110278)
demean = lambda x: x - x.mean()

# De-meaned channel functions
channelFuncsDemean = channelFuncsSubsetPD.transform(demean,axis='columns')

# Plot using lmPlot routine - note this requires conversion to Xarray data type first.
daPlot, daPlotpd, legendList, gFig =  ep.lmPlot(channelFuncsDemean.to_xarray().to_array('t'),
                                                xDim='t', cmap=cmap, mDimLabel='m', thres=1e-4)
Set dataType (No dataType)
Plotting data (No filename), pType=a, thres=0.001, with Seaborn
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
Set dataType (No dataType)
Plotting data (No filename), pType=a, thres=0.0001, with Seaborn
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
In [ ]:
 

FITTING - quick pull from QM3 case-studies¶

Mainly from github/Quantum-Metrology-with-Photoelectrons-Vol3/doc-source/scripts/setup_fit_case-studies_270723.py

NOTE: some things should be cleaned-up and checked.

  1. Is basis not set automatically now...?
  2. Check phases and params more carefully.

Load existing fit data or run fits¶

Note that running fits may be quite time-consuming and computationally intensive, depending on the size of the size of the problem. The default case here will run a small batch for testing if there is no existing data found on the dataPath, otherwise the data is loaded for analysis.

In [19]:
# Configure settings for case study

# Set case study by name
fitSystem='N2O'
#fitStem=f"fit_withNoise_orb5"

# Add noise?
addNoise = 'y'
mu, sigma = 0, 0.05  # Up to approx 10% noise (+/- 0.05)

# Batching - number of fits to run between data dumps
batchSize = 100

# Total fits to run
nMax = 1000

fitStem = f"fit_{orb}_{Eke}"
if addNoise == 'y':
    fitStem += "_noise"
In [20]:
fitStem
Out[20]:
'fit_orb9_8.1_noise'
In [22]:
import xarray as xr
import numpy as np

data.setData('sim', BetaNormX)  # Set simulated data to master structure as "sim"
data.setSubset('sim','AFBLM')   # Set to 'subset' to use for fitting.

# Set basis functions
data.basis = basis

print('\n*Setting  up fit parameters (with constraints)...')

# With auto setting (from existing matrix elements)
data.setMatEFit()

# Report parameter sizes (currently not in main routines)
magCount = 0
phaseCount = 0

for item in data.params:
    if data.params[item].vary:
        if item.startswith('m_'):
            magCount += 1
        else:
            phaseCount +=1
    
print(f"Basis set size = {len(data.params)} params. Fitting with {magCount} magnitudes and {phaseCount} phases floated.")


# Add noise
if addNoise == 'y':
    print(f'\n*** Adding Gaussian noise, mu={mu}, sigma={sigma}')
#     mu, sigma = 0, 0.05  # Up to approx 10% noise (+/- 0.05)
    # creating a noise with the same dimension as the dataset (2,2) 
    noise = np.random.normal(mu, sigma, [data.data['subset']['AFBLM'].t.size, data.data['subset']['AFBLM'].l.size])
    # data.BLMfitPlot()

    # Set noise in Xarray & scale by l
    noiseXR = xr.ones_like(data.data['subset']['AFBLM']) * noise * data.data['subset']['AFBLM'].max()  # FOR OCS ADDED * data.data['subset']['AFBLM'].max() to rescale noise, otherwise ~100%!  Issue is renorm (or not) of values with ADMs.
    # data.data['subset']['AFBLM']['noise'] = ((data.data['subset']['AFBLM'].t, data.data['subset']['AFBLM'].l), noise)
    # xr.where(noiseXR.l>0, noiseXR/noiseXR.l, noiseXR)
    noiseXR = noiseXR.where(noiseXR.l<1, noiseXR/(noiseXR.l))  # Scale by L

    data.data['subset']['AFBLM'] = data.data['subset']['AFBLM'] + noiseXR
    data.data['subset']['AFBLM'] = data.data['subset']['AFBLM'].where(data.data['subset']['AFBLM'].m == 0, 0)
            

print('\n\n*** Setup demo fitting workspace OK.')
Subselected from dataset 'sim', dataType 'AFBLM': 451 from 451 points (100.00%)

*Setting  up fit parameters (with constraints)...
Set 19 complex matrix elements to 38 fitting params, see self.params for details.
Auto-setting parameters.
Parameters
namevalueinitial valueminmaxvaryexpression
m_P_S_P_1_n1_1_1 0.552052840.5520528381595957 1.0000e-04 5.00000000True
m_P_S_P_1_1_n1_1 0.552052840.5520528381595957 1.0000e-04 5.00000000Falsem_P_S_P_1_n1_1_1
m_P_S_P_2_n1_1_1 1.179614001.1796139969198085 1.0000e-04 5.00000000True
m_P_S_P_2_1_n1_1 1.179614001.1796139969198085 1.0000e-04 5.00000000Falsem_P_S_P_2_n1_1_1
m_P_S_P_3_n1_1_1 0.773418890.7734188882567977 1.0000e-04 5.00000000True
m_P_S_P_3_1_n1_1 0.773418890.7734188882567977 1.0000e-04 5.00000000Falsem_P_S_P_3_n1_1_1
m_P_S_P_4_n1_1_1 0.743627800.7436278043671134 1.0000e-04 5.00000000True
m_P_S_P_4_1_n1_1 0.743627800.7436278043671134 1.0000e-04 5.00000000Falsem_P_S_P_4_n1_1_1
m_P_S_P_5_n1_1_1 0.199495560.19949556310846114 1.0000e-04 5.00000000True
m_P_S_P_5_1_n1_1 0.199495560.19949556310846114 1.0000e-04 5.00000000Falsem_P_S_P_5_n1_1_1
m_P_S_P_6_n1_1_1 0.022218890.02221888536181435 1.0000e-04 5.00000000True
m_P_S_P_6_1_n1_1 0.022218890.02221888536181435 1.0000e-04 5.00000000Falsem_P_S_P_6_n1_1_1
m_S_S_S_0_0_0_1 1.829333221.829333218919298 1.0000e-04 5.00000000True
m_S_S_S_1_0_0_1 0.485075870.4850758656806526 1.0000e-04 5.00000000True
m_S_S_S_2_0_0_1 2.489087912.4890879122773666 1.0000e-04 5.00000000True
m_S_S_S_3_0_0_1 0.577691280.5776912775488837 1.0000e-04 5.00000000True
m_S_S_S_4_0_0_1 1.000425721.0004257202501872 1.0000e-04 5.00000000True
m_S_S_S_5_0_0_1 0.236575510.23657550773063152 1.0000e-04 5.00000000True
m_S_S_S_6_0_0_1 0.050528060.05052805754693502 1.0000e-04 5.00000000True
p_P_S_P_1_n1_1_1-2.85739317-2.857393167231562-3.14159265 3.14159265False
p_P_S_P_1_1_n1_1-2.85739317-2.857393167231562-3.14159265 3.14159265Falsep_P_S_P_1_n1_1_1
p_P_S_P_2_n1_1_1 1.843648841.8436488398846065-3.14159265 3.14159265True
p_P_S_P_2_1_n1_1 1.843648841.8436488398846065-3.14159265 3.14159265Falsep_P_S_P_2_n1_1_1
p_P_S_P_3_n1_1_1-0.50485232-0.504852323798376-3.14159265 3.14159265True
p_P_S_P_3_1_n1_1-0.50485232-0.504852323798376-3.14159265 3.14159265Falsep_P_S_P_3_n1_1_1
p_P_S_P_4_n1_1_1-1.39461880-1.3946187996003003-3.14159265 3.14159265True
p_P_S_P_4_1_n1_1-1.39461880-1.3946187996003003-3.14159265 3.14159265Falsep_P_S_P_4_n1_1_1
p_P_S_P_5_n1_1_1 1.889341681.8893416751232721-3.14159265 3.14159265True
p_P_S_P_5_1_n1_1 1.889341681.8893416751232721-3.14159265 3.14159265Falsep_P_S_P_5_n1_1_1
p_P_S_P_6_n1_1_1 2.309593442.309593435686443-3.14159265 3.14159265True
p_P_S_P_6_1_n1_1 2.309593442.309593435686443-3.14159265 3.14159265Falsep_P_S_P_6_n1_1_1
p_S_S_S_0_0_0_1-0.38112270-0.3811227043484284-3.14159265 3.14159265True
p_S_S_S_1_0_0_1 1.609179291.6091792945780479-3.14159265 3.14159265True
p_S_S_S_2_0_0_1-2.17767266-2.177672662861201-3.14159265 3.14159265True
p_S_S_S_3_0_0_1 1.497285401.4972854005089813-3.14159265 3.14159265True
p_S_S_S_4_0_0_1-0.23549828-0.2354982766159596-3.14159265 3.14159265True
p_S_S_S_5_0_0_1 2.108398382.108398381107909-3.14159265 3.14159265True
p_S_S_S_6_0_0_1-2.40162756-2.4016275574535597-3.14159265 3.14159265True
Basis set size = 38 params. Fitting with 13 magnitudes and 12 phases floated.

*** Adding Gaussian noise, mu=0, sigma=0.05


*** Setup demo fitting workspace OK.
In [23]:
# Check cores
# Not sure if using more cores than batch does anything...?
import multiprocessing as mp
nCores = mp.cpu_count()
print(f"nCores={nCores}")

# batchSize = 40
nCores = round(mp.cpu_count() * 0.8) 
num_workers = nCores if (batchSize+2) > nCores else (batchSize+2)

print(f"Setting num_workers={num_workers}")
nCores=56
Setting num_workers=42
In [ ]:
import numpy as np

from datetime import datetime as dt
timeString = dt.now()
print(f"\n*** Running: {timeString.strftime('%Y-%m-%d %H:%M:%S')}")

# Set datapath, 
dataName = 'N2O-fitTest'
dataPath = Path(Path.cwd(),dataName)

# Look for existing Pickle files on path
dataFiles = list(dataPath.expanduser().glob('*.pickle'))

if not dataFiles:
    print("No data found, executing minimal fitting run...")
    
    # Run fit batch - single
    # data.multiFit(nRange = [n,n+batchSize-1], num_workers=batchSize)

    # Run fit batches with checkpoint files
    for n in np.arange(0,nMax,batchSize):
        print(f'*** Running batch [{n},{n+batchSize-1}], {dt.now().strftime("%d%m%y_%H-%M-%S")}')

        # Run fit batch
        data.multiFit(nRange = [n,n+batchSize-1], num_workers=num_workers)

        # Dump data so far - NOTE should pass dir name here!
        data.writeFitData(outStem=f"{fitSystem}_{n+batchSize-1}_{fitStem}")
        
        print(f'Finished batch [{n},{n+batchSize-1}], {dt.now().strftime("%d%m%y_%H-%M-%S")}')
        print(f'Written to file {fitSystem}_{n+batchSize-1}_{fitStem}')

else:
    dataFileIn = dataFiles[-1]   # Add index to select file, although loadFitData will concat multiple files
                                    # Note that concat currently only works for fixed batch sizes however.
    print(f"Set dataFiles: {dataFileIn}")
    data.loadFitData(fList=dataFileIn, dataPath=dataPath)   #.expanduser())
    
    data.BLMfitPlot(keys=['subset','sim'])
    

Post-processing and data overview¶

Post-processing involves aggregation of all the fit run results into a single data structure. This can then be analysed statistically and examined for for best-fit results. In the statistical sense, this is essentailly a search for candidate {{ RADMATE }}, based on the assumption that some of the minima found in the $chi^2$ hyperspace will be the true results. Even if a clear global minima does not exist, searching for candidate {{ RADMATE }} sets based on clustering of results and multiple local minima is still expected to lead to viable candidates provided that the information content of the dataset is sufficient. However, as discussed elsewhere (see {numref}Sect. %s <sect:numerics:fitting-strategies>), in some cases this may not be the case, and other limitations may apply (e.g. certain parameters may be undefined), or additional data required for unique determination of the {{ RADMATE }}.

For more details on the analysis routines, see the {{ PEMtk_docs }}, particularly the fit fidelity and analysis page, and molecular frame analysis data processing page (full analysis for Ref. {cite}hockett2023TopicalReviewExtracting, illustrating the $N_2$ case).

In [25]:
# General stats & post-processing to data tables
data.analyseFits()
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
In [26]:
# The BLMsetPlot routine will output aggregate fit results.
# Here the spread can be taken as a general indication of the uncertainty of 
# the fitting, and indicate whether the fit is well-characterised/the information 
# content of the data is sufficient.
data.BLMsetPlot(xDim = 't', thres=1e-6)  # With xDim and thres set, for more control over outputs
In [36]:
# Write aggregate datasets to HDF5 format
# This is more robust than Pickled data, but PEMtk currently only support output for aggregate (post-processed) fit data.

# 11/04/24 - update dataPath here for current case (manual/slurm datasets)
dataPath = Path(Path.cwd(),dataName + '_1000_080424')

if 'dataFileIn' not in locals().keys():
    if not dataFiles:
        dataFiles = list(dataPath.expanduser().glob('*.pickle'))
        print(dataFiles)
    dataFileIn = dataFiles[-1]

data.processedToHDF5(dataPath = dataPath, outStem = dataFileIn.name, timeStamp=False)
In [45]:
dataFileIn
Out[45]:
PosixPath('/home/jovyan/N2O/epsman2024/proc-fock/N2O-fitTest_1000_080424/N2O_559_fit_orb9_8.1_noise_080424_09-45-35.pickle')
In [55]:
# Check stats...
print(f"Basis set size = {len(data.params)} params. Fitting with {magCount} magnitudes and {phaseCount} phases floated.")
print(f"Data size: {data.data[data.subKey]['ADM'].t.size} t-points")
print(f"Data size: {data.data[data.subKey]['ADM'].size} total ADMs (max (K,Q,S): {data.data[data.subKey]['ADM'].ADM.max().values}).")
print(f"Data size: {data.data[data.subKey]['ADM'].size} total BLMs (max (L,M): {data.data[data.subKey]['AFBLM'].BLM.max().values}).")

if addNoise == 'y':
    print(f"Noise added,  mu={mu}, sigma={sigma}")
Basis set size = 38 params. Fitting with 13 magnitudes and 12 phases floated.
Data size: 41 t-points
Data size: 205 total ADMs (max (K,Q,S): (8, 0, 0)).
Data size: 205 total BLMs (max (L,M): (10, 0)).
Noise added,  mu=0, sigma=0.05
In [73]:
# Histogram fit results (reduced chi^2 vs. fit index)
# This may be quite slow for large datasets, setting limited ranges may help

# Use default auto binning
# data.fitHist()

# Example with range set
# data.fitHist(thres=1.15e-6, bins=100)
data.fitHist(binRange=[1.1004e-6, 1.10045e-6], bins=100)  # TODO: options for dp output in plot?
Mask selected 956 results (from 974).

Here bands in the $\chi^2$ dimension can indicate groupings (local minima) are consistently found. Assuming each grouping is a viable fit candidate parameter set, these can then be explored in further detail.

Test fits:

  • '/home/jovyan/N2O/epsman2024/proc-fock/N2O-fitTest_1000_080424/N2O_559_fit_orb9_8.1_noise_080424_09-45-35.pickle' (1000 fits with noise, 41 t=points), basically see no spread here, all good fits at 1.1e-6.

Data exploration¶

The general aim in this procedure is to ascertain whether there was a good spread of parameters explored, and a single (or few sets) of best-fit results. There are a few procedures and helper methods for this...

View results¶

Single results sets can be viewed in the main data structure, indexed by #.

In [74]:
# Check keys
fitNumber = 696
data.data[fitNumber].keys()
Out[74]:
dict_keys(['AFBLM', 'residual', 'results'])

Here results is an lmFit object, which includes final fit results and information, and AFBLM contains the model (fit) output.

An example is shown below. Of particular note here is which parameters have vary=True - these are included in the fitting - and if there is a column expression, which indicates any parameters defined to have specific relationships (see {numref}Chpt. %s <sect:basis-sets:fitting-intro>). Any correlations found during fitting are also shown, which can also indicate parameters which are related (even if this is not predefined or known a priori).

In [75]:
# Show some results
data.data[fitNumber]['results']
Out[75]:

Fit Result

Fit Statistics
fitting methodleastsq
# function evals2852
# data points451
# variables25
chi-square 4.6876e-04
reduced chi-square 1.1004e-06
Akaike info crit.-6163.37197
Bayesian info crit.-6060.58529
Parameters
namevaluestandard errorrelative errorinitial valueminmaxvaryexpression
m_P_S_P_1_n1_1_1 1.47672449 3.0485e+08(20643663529.67%)0.8860762861939793 1.0000e-04 5.00000000True
m_P_S_P_1_1_n1_1 1.47672449 3.0485e+08(20643663487.72%)0.8860762861939793 1.0000e-04 5.00000000Falsem_P_S_P_1_n1_1_1
m_P_S_P_2_n1_1_1 0.36671056 5.9793e+08(163052950914.54%)0.16069660327581758 1.0000e-04 5.00000000True
m_P_S_P_2_1_n1_1 0.36671056 5.9793e+08(163052950627.77%)0.16069660327581758 1.0000e-04 5.00000000Falsem_P_S_P_2_n1_1_1
m_P_S_P_3_n1_1_1 0.03511796 96545978.0(274919085305.16%)0.620226297126592 1.0000e-04 5.00000000True
m_P_S_P_3_1_n1_1 0.03511796 96545977.8(274919084835.74%)0.620226297126592 1.0000e-04 5.00000000Falsem_P_S_P_3_n1_1_1
m_P_S_P_4_n1_1_1 0.58250076 3.2190e+08(55262131863.93%)0.39760026959116823 1.0000e-04 5.00000000True
m_P_S_P_4_1_n1_1 0.58250076 3.2190e+08(55262131847.24%)0.39760026959116823 1.0000e-04 5.00000000Falsem_P_S_P_4_n1_1_1
m_P_S_P_5_n1_1_1 0.34759066 1.0770e+08(30984570838.56%)0.6241133656185965 1.0000e-04 5.00000000True
m_P_S_P_5_1_n1_1 0.34759066 1.0770e+08(30984570785.71%)0.6241133656185965 1.0000e-04 5.00000000Falsem_P_S_P_5_n1_1_1
m_P_S_P_6_n1_1_1 0.37017292 6.3920e+08(172676967335.07%)0.5467541650735425 1.0000e-04 5.00000000True
m_P_S_P_6_1_n1_1 0.37017292 6.3920e+08(172676967222.45%)0.5467541650735425 1.0000e-04 5.00000000Falsem_P_S_P_6_n1_1_1
m_S_S_S_0_0_0_1 1.52449768 4.7281e+08(31014209552.33%)0.6656859904053511 1.0000e-04 5.00000000True
m_S_S_S_1_0_0_1 1.77896077 1.7710e+08(9955051029.54%)0.5516835821546856 1.0000e-04 5.00000000True
m_S_S_S_2_0_0_1 1.73121799 4.4318e+08(25599061902.89%)0.5518461069601724 1.0000e-04 5.00000000True
m_S_S_S_3_0_0_1 1.14244402 2.8316e+08(24785421198.42%)0.5939359850592106 1.0000e-04 5.00000000True
m_S_S_S_4_0_0_1 1.13376497 2.4828e+08(21898681956.61%)0.17188025192261902 1.0000e-04 5.00000000True
m_S_S_S_5_0_0_1 0.43811343 7.1134e+08(162364836768.44%)0.4455936261077418 1.0000e-04 5.00000000True
m_S_S_S_6_0_0_1 0.16064202 3.7275e+08(232037139803.13%)0.9612133269422948 1.0000e-04 5.00000000True
p_P_S_P_1_n1_1_1-2.85739317 0.00000000(0.00%)-2.857393167231562-3.14159265 3.14159265False
p_P_S_P_1_1_n1_1-2.85739317 0.00000000(0.00%)-2.857393167231562-3.14159265 3.14159265Falsep_P_S_P_1_n1_1_1
p_P_S_P_2_n1_1_1 1.70523025 2.6604e+08(15601541116.80%)0.7700854199580771-3.14159265 3.14159265True
p_P_S_P_2_1_n1_1 1.70523025 2.6604e+08(15601541100.35%)0.7700854199580771-3.14159265 3.14159265Falsep_P_S_P_2_n1_1_1
p_P_S_P_3_n1_1_1-0.23826276 2.6092e+08(109508752144.22%)0.1704705215973239-3.14159265 3.14159265True
p_P_S_P_3_1_n1_1-0.23826276 2.6092e+08(109508752340.72%)0.1704705215973239-3.14159265 3.14159265Falsep_P_S_P_3_n1_1_1
p_P_S_P_4_n1_1_1-0.32836413 3.0618e+08(93244875233.47%)0.5035801054886633-3.14159265 3.14159265True
p_P_S_P_4_1_n1_1-0.32836413 3.0618e+08(93244875017.84%)0.5035801054886633-3.14159265 3.14159265Falsep_P_S_P_4_n1_1_1
p_P_S_P_5_n1_1_1 0.29479429 3.3014e+08(111990060820.57%)0.23470718604773666-3.14159265 3.14159265True
p_P_S_P_5_1_n1_1 0.29479429 3.3014e+08(111990061073.08%)0.23470718604773666-3.14159265 3.14159265Falsep_P_S_P_5_n1_1_1
p_P_S_P_6_n1_1_1 3.12014299 4.9044e+08(15718633629.78%)0.7051031845542469-3.14159265 3.14159265True
p_P_S_P_6_1_n1_1 3.12014299 4.9044e+08(15718633605.37%)0.7051031845542469-3.14159265 3.14159265Falsep_P_S_P_6_n1_1_1
p_S_S_S_0_0_0_1 1.28842864 5.9702e+08(46337338605.79%)0.43063262562515603-3.14159265 3.14159265True
p_S_S_S_1_0_0_1 0.28338855 1.7468e+08(61638350418.97%)0.7355325030452876-3.14159265 3.14159265True
p_S_S_S_2_0_0_1-1.46879940 4.6032e+08(31339918828.80%)0.06889941894672202-3.14159265 3.14159265True
p_S_S_S_3_0_0_1-0.17951447 1.2102e+08(67416928787.13%)0.033275669140320874-3.14159265 3.14159265True
p_S_S_S_4_0_0_1 1.02789099 5.0388e+08(49020502367.99%)0.3138781475801862-3.14159265 3.14159265True
p_S_S_S_5_0_0_1 1.47848227 4.4368e+08(30008931072.01%)0.7675345196637045-3.14159265 3.14159265True
p_S_S_S_6_0_0_1-0.30925984 3.4011e+08(109976955515.33%)0.6836014802961643-3.14159265 3.14159265True
Correlations (unreported values are < 0.100)
Parameter1Parameter 2Correlation
m_P_S_P_4_n1_1_1m_P_S_P_6_n1_1_1+0.9500
m_S_S_S_5_0_0_1m_S_S_S_6_0_0_1-0.9096
m_P_S_P_1_n1_1_1m_P_S_P_6_n1_1_1-0.9071
m_S_S_S_3_0_0_1m_S_S_S_4_0_0_1+0.8774
m_P_S_P_1_n1_1_1p_S_S_S_2_0_0_1+0.8649
m_P_S_P_1_n1_1_1m_P_S_P_4_n1_1_1-0.8594
m_P_S_P_6_n1_1_1p_S_S_S_2_0_0_1-0.8498
m_P_S_P_2_n1_1_1m_P_S_P_3_n1_1_1+0.8411
m_P_S_P_5_n1_1_1m_S_S_S_3_0_0_1+0.8211
m_P_S_P_4_n1_1_1p_S_S_S_2_0_0_1-0.8113
m_P_S_P_5_n1_1_1m_P_S_P_6_n1_1_1+0.7894
m_S_S_S_5_0_0_1p_S_S_S_4_0_0_1+0.7710
p_P_S_P_2_n1_1_1p_S_S_S_2_0_0_1+0.7649
m_P_S_P_5_n1_1_1m_S_S_S_4_0_0_1+0.7610
m_P_S_P_1_n1_1_1m_P_S_P_3_n1_1_1-0.7121
m_S_S_S_1_0_0_1p_P_S_P_6_n1_1_1+0.7042
m_P_S_P_4_n1_1_1m_P_S_P_5_n1_1_1+0.7038
m_P_S_P_4_n1_1_1m_S_S_S_3_0_0_1+0.7006
m_S_S_S_6_0_0_1p_S_S_S_4_0_0_1-0.7005
m_P_S_P_6_n1_1_1m_S_S_S_3_0_0_1+0.6814
m_P_S_P_5_n1_1_1p_S_S_S_4_0_0_1-0.6652
m_P_S_P_2_n1_1_1m_P_S_P_5_n1_1_1-0.6604
m_S_S_S_1_0_0_1p_S_S_S_5_0_0_1-0.6575
m_P_S_P_6_n1_1_1p_S_S_S_4_0_0_1-0.6551
m_P_S_P_2_n1_1_1m_S_S_S_3_0_0_1-0.6378
m_P_S_P_3_n1_1_1p_P_S_P_2_n1_1_1-0.6287
m_S_S_S_4_0_0_1p_S_S_S_0_0_0_1+0.5975
m_P_S_P_6_n1_1_1m_S_S_S_4_0_0_1+0.5922
p_P_S_P_6_n1_1_1p_S_S_S_5_0_0_1-0.5856
m_S_S_S_3_0_0_1p_S_S_S_2_0_0_1-0.5844
m_P_S_P_2_n1_1_1p_S_S_S_0_0_0_1-0.5721
m_S_S_S_0_0_0_1m_S_S_S_6_0_0_1+0.5679
m_S_S_S_0_0_0_1p_S_S_S_5_0_0_1+0.5630
m_P_S_P_4_n1_1_1p_S_S_S_4_0_0_1-0.5612
m_S_S_S_4_0_0_1p_S_S_S_2_0_0_1-0.5611
m_P_S_P_5_n1_1_1p_S_S_S_2_0_0_1-0.5585
m_P_S_P_3_n1_1_1m_S_S_S_1_0_0_1+0.5500
m_P_S_P_4_n1_1_1m_S_S_S_4_0_0_1+0.5495
m_S_S_S_0_0_0_1m_S_S_S_5_0_0_1-0.5479
m_P_S_P_5_n1_1_1p_S_S_S_0_0_0_1+0.5460
m_P_S_P_3_n1_1_1p_S_S_S_2_0_0_1-0.5404
m_S_S_S_0_0_0_1m_S_S_S_2_0_0_1-0.5336
m_S_S_S_3_0_0_1p_S_S_S_0_0_0_1+0.5201
m_P_S_P_1_n1_1_1p_P_S_P_2_n1_1_1+0.5138
m_P_S_P_2_n1_1_1m_S_S_S_4_0_0_1-0.5108
m_S_S_S_6_0_0_1p_S_S_S_0_0_0_1-0.5083
m_P_S_P_1_n1_1_1p_S_S_S_4_0_0_1+0.5077
m_S_S_S_2_0_0_1m_S_S_S_4_0_0_1-0.5068
m_P_S_P_2_n1_1_1m_S_S_S_1_0_0_1+0.5065
m_S_S_S_0_0_0_1p_S_S_S_1_0_0_1-0.5003
m_S_S_S_1_0_0_1p_S_S_S_6_0_0_1-0.4942
p_S_S_S_2_0_0_1p_S_S_S_4_0_0_1+0.4902
m_S_S_S_3_0_0_1p_S_S_S_3_0_0_1-0.4901
p_S_S_S_1_0_0_1p_S_S_S_3_0_0_1-0.4846
m_P_S_P_4_n1_1_1p_S_S_S_0_0_0_1+0.4838
m_S_S_S_0_0_0_1p_S_S_S_0_0_0_1-0.4736
m_S_S_S_1_0_0_1m_S_S_S_2_0_0_1-0.4725
m_P_S_P_1_n1_1_1m_P_S_P_5_n1_1_1-0.4685
m_S_S_S_6_0_0_1p_S_S_S_5_0_0_1+0.4638
m_S_S_S_1_0_0_1p_S_S_S_4_0_0_1+0.4584
m_S_S_S_3_0_0_1p_S_S_S_1_0_0_1+0.4455
p_P_S_P_3_n1_1_1p_P_S_P_4_n1_1_1-0.4408
m_P_S_P_5_n1_1_1p_P_S_P_6_n1_1_1-0.4401
p_P_S_P_4_n1_1_1p_S_S_S_5_0_0_1+0.4320
m_P_S_P_4_n1_1_1m_S_S_S_5_0_0_1-0.4313
p_P_S_P_2_n1_1_1p_P_S_P_4_n1_1_1+0.4312
m_P_S_P_3_n1_1_1p_S_S_S_5_0_0_1-0.4310
m_P_S_P_2_n1_1_1p_P_S_P_2_n1_1_1-0.4229
m_P_S_P_6_n1_1_1p_S_S_S_0_0_0_1+0.4093
m_P_S_P_6_n1_1_1m_S_S_S_5_0_0_1-0.4091
m_P_S_P_1_n1_1_1m_S_S_S_3_0_0_1-0.4063
m_P_S_P_1_n1_1_1m_S_S_S_1_0_0_1-0.3961
p_P_S_P_4_n1_1_1p_P_S_P_6_n1_1_1-0.3952
p_P_S_P_4_n1_1_1p_S_S_S_2_0_0_1+0.3909
m_S_S_S_6_0_0_1p_S_S_S_6_0_0_1+0.3906
m_P_S_P_1_n1_1_1m_S_S_S_6_0_0_1-0.3851
m_P_S_P_2_n1_1_1p_P_S_P_6_n1_1_1+0.3833
m_P_S_P_1_n1_1_1p_S_S_S_5_0_0_1+0.3817
p_S_S_S_0_0_0_1p_S_S_S_5_0_0_1-0.3779
m_P_S_P_3_n1_1_1m_P_S_P_6_n1_1_1+0.3753
m_P_S_P_6_n1_1_1m_S_S_S_6_0_0_1+0.3747
m_S_S_S_3_0_0_1p_S_S_S_4_0_0_1-0.3738
p_P_S_P_6_n1_1_1p_S_S_S_6_0_0_1-0.3716
m_P_S_P_1_n1_1_1m_S_S_S_5_0_0_1+0.3648
m_S_S_S_1_0_0_1m_S_S_S_6_0_0_1-0.3556
m_P_S_P_1_n1_1_1m_S_S_S_4_0_0_1-0.3540
m_S_S_S_5_0_0_1p_S_S_S_5_0_0_1-0.3528
p_P_S_P_6_n1_1_1p_S_S_S_4_0_0_1+0.3510
m_P_S_P_6_n1_1_1p_P_S_P_2_n1_1_1-0.3503
p_S_S_S_2_0_0_1p_S_S_S_5_0_0_1+0.3489
p_P_S_P_4_n1_1_1p_P_S_P_5_n1_1_1-0.3418
m_S_S_S_1_0_0_1m_S_S_S_5_0_0_1+0.3406
m_S_S_S_0_0_0_1p_S_S_S_3_0_0_1+0.3355
m_S_S_S_4_0_0_1p_S_S_S_3_0_0_1-0.3349
m_P_S_P_3_n1_1_1p_P_S_P_6_n1_1_1+0.3339
p_S_S_S_5_0_0_1p_S_S_S_6_0_0_1+0.3304
m_P_S_P_2_n1_1_1p_S_S_S_5_0_0_1-0.3247
p_P_S_P_4_n1_1_1p_S_S_S_1_0_0_1-0.3231
m_P_S_P_5_n1_1_1m_S_S_S_1_0_0_1-0.3218
m_S_S_S_5_0_0_1p_S_S_S_0_0_0_1+0.3212
m_P_S_P_4_n1_1_1m_S_S_S_6_0_0_1+0.3175
m_S_S_S_2_0_0_1m_S_S_S_3_0_0_1-0.3163
m_P_S_P_4_n1_1_1p_P_S_P_4_n1_1_1-0.3122
m_P_S_P_4_n1_1_1p_S_S_S_5_0_0_1-0.3120
m_S_S_S_0_0_0_1m_S_S_S_3_0_0_1-0.3116
p_P_S_P_2_n1_1_1p_S_S_S_5_0_0_1+0.3093
m_S_S_S_4_0_0_1p_P_S_P_6_n1_1_1-0.3021
m_P_S_P_4_n1_1_1m_S_S_S_2_0_0_1-0.3012
p_P_S_P_2_n1_1_1p_S_S_S_0_0_0_1+0.3004
m_S_S_S_0_0_0_1m_S_S_S_4_0_0_1-0.2972
m_S_S_S_1_0_0_1p_P_S_P_4_n1_1_1-0.2969
m_P_S_P_1_n1_1_1m_S_S_S_2_0_0_1+0.2889
m_P_S_P_3_n1_1_1p_S_S_S_0_0_0_1-0.2884
m_S_S_S_1_0_0_1p_S_S_S_2_0_0_1-0.2883
m_P_S_P_6_n1_1_1m_S_S_S_2_0_0_1-0.2866
m_P_S_P_5_n1_1_1m_S_S_S_5_0_0_1-0.2827
m_P_S_P_1_n1_1_1p_P_S_P_6_n1_1_1-0.2816
m_P_S_P_3_n1_1_1p_P_S_P_3_n1_1_1-0.2779
m_P_S_P_4_n1_1_1p_P_S_P_2_n1_1_1-0.2727
m_S_S_S_2_0_0_1p_S_S_S_2_0_0_1+0.2725
m_S_S_S_1_0_0_1p_S_S_S_0_0_0_1+0.2705
m_P_S_P_3_n1_1_1m_P_S_P_4_n1_1_1+0.2705
p_P_S_P_2_n1_1_1p_S_S_S_4_0_0_1+0.2674
p_S_S_S_4_0_0_1p_S_S_S_5_0_0_1-0.2642
m_S_S_S_0_0_0_1p_S_S_S_2_0_0_1+0.2638
m_S_S_S_2_0_0_1p_S_S_S_1_0_0_1+0.2636
m_P_S_P_4_n1_1_1p_P_S_P_6_n1_1_1+0.2596
m_P_S_P_2_n1_1_1p_S_S_S_4_0_0_1+0.2592
m_P_S_P_1_n1_1_1m_P_S_P_2_n1_1_1-0.2487
m_S_S_S_6_0_0_1p_P_S_P_5_n1_1_1-0.2482
m_P_S_P_2_n1_1_1m_P_S_P_4_n1_1_1-0.2404
m_P_S_P_5_n1_1_1m_S_S_S_0_0_0_1-0.2387
m_P_S_P_4_n1_1_1m_S_S_S_1_0_0_1+0.2384
m_S_S_S_2_0_0_1p_S_S_S_5_0_0_1-0.2352
m_S_S_S_1_0_0_1p_S_S_S_3_0_0_1+0.2336
m_S_S_S_4_0_0_1p_S_S_S_1_0_0_1+0.2321
p_S_S_S_0_0_0_1p_S_S_S_6_0_0_1-0.2296
m_S_S_S_1_0_0_1p_P_S_P_2_n1_1_1-0.2236
m_S_S_S_3_0_0_1p_P_S_P_4_n1_1_1-0.2214
m_P_S_P_1_n1_1_1p_P_S_P_3_n1_1_1+0.2211
m_P_S_P_5_n1_1_1m_S_S_S_6_0_0_1+0.2203
m_S_S_S_4_0_0_1p_P_S_P_5_n1_1_1+0.2186
m_S_S_S_6_0_0_1p_S_S_S_2_0_0_1-0.2172
m_S_S_S_0_0_0_1p_P_S_P_2_n1_1_1+0.2147
m_P_S_P_5_n1_1_1p_S_S_S_6_0_0_1+0.2081
m_P_S_P_3_n1_1_1m_S_S_S_3_0_0_1-0.2055
m_P_S_P_2_n1_1_1m_S_S_S_0_0_0_1+0.2050
p_P_S_P_4_n1_1_1p_S_S_S_3_0_0_1+0.2047
m_S_S_S_3_0_0_1p_S_S_S_6_0_0_1+0.2041
m_P_S_P_3_n1_1_1m_S_S_S_6_0_0_1+0.2040
m_S_S_S_1_0_0_1p_S_S_S_1_0_0_1-0.2039
p_S_S_S_1_0_0_1p_S_S_S_2_0_0_1-0.2030
m_S_S_S_0_0_0_1p_P_S_P_4_n1_1_1+0.1992
m_P_S_P_6_n1_1_1p_S_S_S_5_0_0_1-0.1992
p_S_S_S_0_0_0_1p_S_S_S_2_0_0_1-0.1984
m_S_S_S_5_0_0_1p_S_S_S_6_0_0_1-0.1971
m_S_S_S_2_0_0_1p_S_S_S_0_0_0_1-0.1956
p_P_S_P_2_n1_1_1p_P_S_P_5_n1_1_1-0.1948
m_S_S_S_3_0_0_1p_P_S_P_6_n1_1_1-0.1930
p_S_S_S_3_0_0_1p_S_S_S_6_0_0_1-0.1918
p_P_S_P_3_n1_1_1p_P_S_P_5_n1_1_1+0.1914
m_S_S_S_1_0_0_1m_S_S_S_3_0_0_1-0.1890
m_S_S_S_4_0_0_1p_S_S_S_4_0_0_1-0.1884
p_S_S_S_4_0_0_1p_S_S_S_6_0_0_1-0.1824
m_P_S_P_6_n1_1_1p_P_S_P_3_n1_1_1-0.1815
m_P_S_P_1_n1_1_1p_S_S_S_0_0_0_1-0.1814
m_S_S_S_4_0_0_1p_P_S_P_2_n1_1_1-0.1793
p_P_S_P_2_n1_1_1p_P_S_P_3_n1_1_1-0.1783
m_P_S_P_3_n1_1_1m_P_S_P_5_n1_1_1-0.1759
m_S_S_S_6_0_0_1p_P_S_P_3_n1_1_1-0.1722
m_P_S_P_4_n1_1_1p_S_S_S_1_0_0_1+0.1715
m_P_S_P_5_n1_1_1m_S_S_S_2_0_0_1-0.1711
m_P_S_P_5_n1_1_1p_S_S_S_5_0_0_1+0.1711
m_P_S_P_2_n1_1_1m_P_S_P_6_n1_1_1-0.1697
m_S_S_S_5_0_0_1p_S_S_S_2_0_0_1+0.1650
m_S_S_S_4_0_0_1m_S_S_S_6_0_0_1-0.1649
m_P_S_P_1_n1_1_1p_P_S_P_5_n1_1_1+0.1644
m_S_S_S_3_0_0_1m_S_S_S_5_0_0_1-0.1641
m_P_S_P_3_n1_1_1p_P_S_P_5_n1_1_1-0.1638
m_P_S_P_2_n1_1_1p_P_S_P_5_n1_1_1-0.1618
m_S_S_S_2_0_0_1p_P_S_P_6_n1_1_1-0.1615
m_S_S_S_0_0_0_1p_S_S_S_4_0_0_1-0.1613
m_P_S_P_2_n1_1_1p_P_S_P_3_n1_1_1-0.1567
p_S_S_S_2_0_0_1p_S_S_S_3_0_0_1+0.1545
m_P_S_P_2_n1_1_1p_S_S_S_1_0_0_1-0.1543
m_P_S_P_4_n1_1_1m_S_S_S_0_0_0_1-0.1521
m_P_S_P_2_n1_1_1p_S_S_S_6_0_0_1-0.1512
m_S_S_S_4_0_0_1m_S_S_S_5_0_0_1+0.1510
m_P_S_P_6_n1_1_1m_S_S_S_0_0_0_1-0.1466
m_P_S_P_3_n1_1_1m_S_S_S_2_0_0_1-0.1466
m_S_S_S_1_0_0_1p_P_S_P_5_n1_1_1-0.1401
m_P_S_P_6_n1_1_1m_S_S_S_1_0_0_1+0.1388
p_P_S_P_5_n1_1_1p_S_S_S_3_0_0_1-0.1386
p_P_S_P_2_n1_1_1p_S_S_S_1_0_0_1-0.1351
p_P_S_P_3_n1_1_1p_S_S_S_6_0_0_1-0.1323
p_P_S_P_2_n1_1_1p_P_S_P_6_n1_1_1+0.1296
m_S_S_S_3_0_0_1m_S_S_S_6_0_0_1+0.1288
m_P_S_P_4_n1_1_1p_S_S_S_3_0_0_1-0.1233
m_P_S_P_6_n1_1_1p_P_S_P_5_n1_1_1-0.1218
p_P_S_P_2_n1_1_1p_S_S_S_6_0_0_1-0.1205
p_S_S_S_0_0_0_1p_S_S_S_4_0_0_1+0.1199
m_S_S_S_3_0_0_1p_S_S_S_5_0_0_1+0.1189
m_P_S_P_1_n1_1_1p_P_S_P_4_n1_1_1+0.1174
m_S_S_S_2_0_0_1p_S_S_S_6_0_0_1+0.1173
m_S_S_S_0_0_0_1p_P_S_P_6_n1_1_1+0.1152
m_S_S_S_3_0_0_1p_P_S_P_2_n1_1_1-0.1117
m_P_S_P_2_n1_1_1m_S_S_S_5_0_0_1+0.1102
p_P_S_P_6_n1_1_1p_S_S_S_0_0_0_1+0.1101
m_P_S_P_2_n1_1_1p_S_S_S_3_0_0_1+0.1079
p_P_S_P_2_n1_1_1p_S_S_S_3_0_0_1+0.1051
m_P_S_P_6_n1_1_1p_S_S_S_1_0_0_1+0.1025
m_P_S_P_5_n1_1_1p_P_S_P_4_n1_1_1+0.1021
p_P_S_P_4_n1_1_1p_S_S_S_6_0_0_1+0.1019
m_P_S_P_5_n1_1_1p_S_S_S_1_0_0_1+0.1019
m_P_S_P_3_n1_1_1m_S_S_S_4_0_0_1-0.1017
m_S_S_S_5_0_0_1p_P_S_P_6_n1_1_1-0.1002

Classify candidate sets¶

To probe the minima found, the classifyFits method can be used. This bins results into "candidate" groups, which can then be examined in detail.

In [86]:
# Run with defaults
# data.classifyFits()

# For more control, pass bins
# Here the minima is set at one end, and a %age range used for bins
minVal = data.fitsSummary['Stats']['redchi']['min']    

# binRangePC = 0.9e-6   # Gives ~3 per bin for first few best fits
binRangePC = 1.5e-6   # Gives ~6 per bin

data.classifyFits(bins = [minVal, minVal + binRangePC*minVal , 20])
success chisqr redchi
count unique top freq count unique top freq count unique top freq
redchiGroup
A 6 1 True 6 6.0 6.0 0.000469 1.0 6.0 6.0 0.000001 1.0
B 4 1 True 4 4.0 4.0 0.000469 1.0 4.0 4.0 0.000001 1.0
C 6 1 True 6 6.0 6.0 0.000469 1.0 6.0 6.0 0.000001 1.0
D 8 1 True 8 8.0 8.0 0.000469 1.0 8.0 8.0 0.000001 1.0
E 10 1 True 10 10.0 10.0 0.000469 1.0 10.0 10.0 0.000001 1.0
F 10 1 True 10 10.0 10.0 0.000469 1.0 10.0 10.0 0.000001 1.0
G 15 1 True 15 15.0 15.0 0.000469 1.0 15.0 15.0 0.000001 1.0
H 16 1 True 16 16.0 16.0 0.000469 1.0 16.0 16.0 0.000001 1.0
I 15 1 True 15 15.0 15.0 0.000469 1.0 15.0 15.0 0.000001 1.0
J 15 1 True 15 15.0 15.0 0.000469 1.0 15.0 15.0 0.000001 1.0
K 25 1 True 25 25.0 25.0 0.000469 1.0 25.0 25.0 0.000001 1.0
L 17 1 True 17 17.0 17.0 0.000469 1.0 17.0 17.0 0.000001 1.0
M 23 1 True 23 23.0 23.0 0.000469 1.0 23.0 23.0 0.000001 1.0
N 28 1 True 28 28.0 28.0 0.000469 1.0 28.0 28.0 0.000001 1.0
O 17 1 True 17 17.0 17.0 0.000469 1.0 17.0 17.0 0.000001 1.0
P 19 1 True 19 19.0 19.0 0.000469 1.0 19.0 19.0 0.000001 1.0
Q 11 1 True 11 11.0 11.0 0.000469 1.0 11.0 11.0 0.000001 1.0
R 22 1 True 22 22.0 22.0 0.000469 1.0 22.0 22.0 0.000001 1.0
S 14 1 True 14 14.0 14.0 0.000469 1.0 14.0 14.0 0.000001 1.0
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

Explore candidate result sets¶

Drill-down on a candidate set of results, and examine values and spreads. For more details see {{ PEMtk_docs }}, especially the analysis routines page. (See also {numref}Sect. %s <sect:platform:pythonEcosystem> for details on the plotting libaries implemented here.)

Raw results¶

Plot spreads in magnitude and phase parameters. Statistical plots are available for Seaborn and Holoviews backends, with some slightly different options.

In [87]:
# From the candidates, select a group for analysis
selGroup = 'A'
In [88]:
# paramPlot can be used to check the spread on each parameter.
# Plots use Seaborn or Holoviews/Bokeh
# Colour-mapping is controlled by the 'hue' paramter, additionally pass hRound for sig. fig control.
# The remap setting allows for short-hand labels as set in data.lmmu

paramType = 'm' # Set for (m)agnitude or (p)hase parameters
hRound = 14 # Set for cmapping, default may be too small (leads to all grey cmap on points)

paramPlotBackend = 'hv'  # 'hv' or 'sns'

data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', 
               backend=paramPlotBackend, hvType='violin', 
               returnFlag = True, hRound=hRound, remap = 'lmMap');
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
In [89]:
paramType = 'p' # Set for (m)agnitude or (p)hase parameters
data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', backend=paramPlotBackend, hvType='violin', 
               returnFlag = True, hRound=hRound, remap = 'lmMap'); 
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

Phases, phase shifts & corrections¶

Depending on how the fit was configured, phases may be defined in different ways. To set the phases relative to a speific parameter, and wrap to a specified range, use the phaseCorrection() method. This defaults to using the first parameter as a reference phase, and wraps to $-\pi:\pi$. The phase-corrected values are output to a new Type, 'pc', and a set of normalised magnitudes to 'n'. Additional settings can be passed for more control, as shown below.

In [90]:
# Run phase correction routine
# Set absFlag=True for unsigned phases (mapped to 0:pi)
# Set useRef=False to set ref phase as 0, otherwise the reference value is set.
phaseCorrParams={'absFlag':True, 'useRef':False}
data.phaseCorrection(**phaseCorrParams)  
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
Set ref param = P_S_P_1_1_n1_1
Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 P_S_P_6_1_n1_1 P_S_P_6_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1
Fit Type redchiGroup
0 m O 0.587109 0.587109 0.876792 0.876792 0.826907 0.826907 0.819588 0.819588 0.557993 0.557993 0.364997 0.364997 1.749348 1.110452 1.540857 1.284653 1.464809 0.892935 0.182733
n O 0.141889 0.141889 0.211898 0.211898 0.199842 0.199842 0.198074 0.198074 0.134853 0.134853 0.088211 0.088211 0.422773 0.268368 0.372386 0.310468 0.354007 0.215800 0.044162
p O -2.857393 -2.857393 1.473154 1.473154 2.109622 2.109622 -1.844451 -1.844451 -0.706675 -0.706675 1.182573 1.182573 2.511460 0.363996 0.229466 0.475833 2.130921 2.365797 -1.358567
pc O 0.000000 0.000000 1.952638 1.952638 1.316170 1.316170 1.012942 1.012942 2.150718 2.150718 2.243219 2.243219 0.914332 3.061796 3.086860 2.949959 1.294871 1.059995 1.498826
2 m P 1.465974 1.465974 0.374354 0.374354 0.165118 0.165118 0.573602 0.573602 0.356709 0.356709 0.375444 0.375444 2.088999 1.315278 0.440169 1.114193 1.746866 0.825574 0.178978
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
994 pc J 0.000000 0.000000 1.801854 1.801854 2.748622 2.748622 0.307106 0.307106 3.106621 3.106621 2.966311 2.966311 1.274185 2.988401 1.494340 2.859567 1.031810 1.642137 0.776052
997 m A 1.071010 1.071010 1.058500 1.058500 0.217025 0.217025 0.665377 0.665377 0.147383 0.147383 0.366794 0.366794 1.085956 1.453644 2.060101 1.511633 1.135005 0.398858 0.159724
n A 0.258835 0.258835 0.255812 0.255812 0.052449 0.052449 0.160804 0.160804 0.035619 0.035619 0.088645 0.088645 0.262448 0.351308 0.497873 0.365323 0.274302 0.096394 0.038601
p A -2.857393 -2.857393 -1.464188 -1.464188 0.428979 0.428979 0.890176 0.890176 -0.526859 -0.526859 3.141593 3.141593 -0.945850 -0.661391 2.185180 -1.253769 -1.084264 -0.028552 0.198967
pc A 0.000000 0.000000 1.393205 1.393205 2.996813 2.996813 2.535616 2.535616 2.330534 2.330534 0.284199 0.284199 1.911543 2.196002 1.240612 1.603625 1.773129 2.828841 3.056360

1124 rows × 19 columns

Examine new data types...

In [91]:
paramType = 'n'
data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', 
               backend=paramPlotBackend, hvType='violin', kind='box',
               returnFlag = True, hRound=hRound, remap = 'lmMap');
WARNING:param.Scatter10901: Setting non-parameter attribute kind=box using a mechanism intended only for parameters
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
In [92]:
paramType = 'pc'
data.paramPlot(selectors={'Type':paramType, 'redchiGroup':selGroup}, hue = 'redchi', 
               backend=paramPlotBackend, hvType='violin', kind='box',
               returnFlag = True, hRound=hRound, remap = 'lmMap');
WARNING:param.Scatter11608: Setting non-parameter attribute kind=box using a mechanism intended only for parameters
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.

Parameter estimation & fidelity¶

For case studies, the fit results can be directly compared to the known input parameters. This should give a feel for how well the data defines the matrix elements (parameters) in this case. In general, probing the correlations and spread of results, and comparing to other (unfitted) results is required to estimate fidelity, see {{ QM12 }} for further discussion.

Best values and statistics¶

To get a final parameter set and associated statistics, based on a subset of the fit results, the paramsReport() method is available. If reference data is available, as for the case studies herein, the paramsCompare() method can also be used to compare with the reference case.

In [93]:
# Parameter summary
data.paramsReport(inds = {'redchiGroup':selGroup})
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
In [94]:
# Parameter comparison
# Note this uses phaseCorrParams as set previously for consistency
data.paramsCompare(phaseCorrParams=phaseCorrParams)
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
Set ref param = P_S_P_1_1_n1_1
Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 P_S_P_6_1_n1_1 P_S_P_6_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1
Fit Type
ref m 0.552053 0.552053 1.179614 1.179614 0.773419 0.773419 0.743628 0.743628 0.199496 0.199496 0.022219 0.022219 1.829333 0.485076 2.489088 0.577691 1.000426 0.236576 0.050528
n 0.134106 0.134106 0.286555 0.286555 0.187881 0.187881 0.180644 0.180644 0.048462 0.048462 0.005397 0.005397 0.444387 0.117836 0.604657 0.140335 0.243026 0.057470 0.012274
p -2.857393 -2.857393 1.843649 1.843649 -0.504852 -0.504852 -1.394619 -1.394619 1.889342 1.889342 2.309593 2.309593 -0.381123 1.609179 -2.177673 1.497285 -0.235498 2.108398 -2.401628
pc 0.000000 0.000000 1.582143 1.582143 2.352541 2.352541 1.462774 1.462774 1.536450 1.536450 1.116199 1.116199 2.476270 1.816613 0.679721 1.928507 2.621895 1.317394 0.455766
In [95]:
# Display above results With column name remapping to (l,m) labels only

# With Pandas functionality
data.paramsSummaryComp.rename(columns=data.lmmu['lmMap'])

# With utility method
# summaryRenamed = pemtk.fit._util.renameParams(data.paramsSummaryComp, data.lmmu['lmMap']) 
# summaryRenamed
Out[95]:
Param 1,1 1,-1 2,1 2,-1 3,1 3,-1 4,1 4,-1 5,1 5,-1 6,1 6,-1 0,0 1,0 2,0 3,0 4,0 5,0 6,0
Type Source dType
m mean num 1.055450 1.055450 0.754709 0.754709 0.444572 0.444572 0.697056 0.697056 0.384742 0.384742 0.365072 0.365072 1.650705 1.434879 1.513031 1.356301 1.345129 0.483148 0.150956
ref num 0.552053 0.552053 1.179614 1.179614 0.773419 0.773419 0.743628 0.743628 0.199496 0.199496 0.022219 0.022219 1.829333 0.485076 2.489088 0.577691 1.000426 0.236576 0.050528
diff % 47.695018 47.695018 56.300416 56.300416 73.969488 73.969488 6.681234 6.681234 48.148163 48.148163 93.913830 93.913830 10.821361 66.193949 64.510079 57.406851 25.626034 51.034553 66.527989
num 0.503397 0.503397 -0.424905 -0.424905 -0.328847 -0.328847 -0.046572 -0.046572 0.185246 0.185246 0.342853 0.342853 -0.178629 0.949803 -0.976057 0.778610 0.344703 0.246572 0.100428
std % 29.946654 29.946654 45.293117 45.293117 62.321106 62.321106 12.117068 12.117068 34.344714 34.344714 2.168152 2.168152 21.148320 14.386940 22.936991 10.259865 14.721248 39.056815 11.613643
num 0.316072 0.316072 0.341831 0.341831 0.277062 0.277062 0.084463 0.084463 0.132138 0.132138 0.007915 0.007915 0.349096 0.206435 0.347044 0.139155 0.198020 0.188702 0.017532
diff/std % 159.266602 159.266602 124.302364 124.302364 118.690911 118.690911 55.139027 55.139027 140.190898 140.190898 4331.515603 4331.515603 51.168893 460.097493 281.249091 559.528328 174.075143 130.667472 572.843408
n mean num 0.255075 0.255075 0.182394 0.182394 0.107442 0.107442 0.168461 0.168461 0.092982 0.092982 0.088229 0.088229 0.398933 0.346774 0.365661 0.327783 0.325083 0.116765 0.036482
ref num 0.134106 0.134106 0.286555 0.286555 0.187881 0.187881 0.180644 0.180644 0.048462 0.048462 0.005397 0.005397 0.444387 0.117836 0.604657 0.140335 0.243026 0.057470 0.012274
diff % 47.424829 47.424829 57.107972 57.107972 74.868228 74.868228 7.232368 7.232368 47.880303 47.880303 93.882388 93.882388 11.393844 66.019301 65.359959 57.186805 25.241811 50.781629 66.355067
num 0.120969 0.120969 -0.104161 -0.104161 -0.080440 -0.080440 -0.012184 -0.012184 0.044520 0.044520 0.082831 0.082831 -0.045454 0.228938 -0.238996 0.187449 0.082057 0.059295 0.024208
std % 29.946769 29.946769 45.292977 45.292977 62.321130 62.321130 12.116934 12.116934 34.344739 34.344739 2.167948 2.167948 21.148450 14.386846 22.936949 10.259669 14.721217 39.057007 11.613524
num 0.076387 0.076387 0.082612 0.082612 0.066959 0.066959 0.020412 0.020412 0.031935 0.031935 0.001913 0.001913 0.084368 0.049890 0.083871 0.033629 0.047856 0.045605 0.004237
diff/std % 158.363758 158.363758 126.085712 126.085712 120.132976 120.132976 59.688102 59.688102 139.410879 139.410879 4330.472681 4330.472681 53.875552 458.886539 284.954904 557.394275 171.465524 130.019254 571.360323
p mean num -2.857393 -2.857393 -0.462299 -0.462299 0.413259 0.413259 1.358097 1.358097 0.292101 0.292101 1.187404 1.187404 0.282125 0.485615 0.400679 0.349399 0.217846 2.145848 -0.913459
ref num -2.857393 -2.857393 1.843649 1.843649 -0.504852 -0.504852 -1.394619 -1.394619 1.889342 1.889342 2.309593 2.309593 -0.381123 1.609179 -2.177673 1.497285 -0.235498 2.108398 -2.401628
diff % 0.000000 0.000000 498.800404 498.800404 222.163669 222.163669 202.689219 202.689219 546.811368 546.811368 94.507767 94.507767 235.090249 231.369023 643.495002 328.532087 208.103055 1.745235 162.915649
num 0.000000 0.000000 -2.305947 -2.305947 0.918111 0.918111 2.752715 2.752715 -1.597241 -1.597241 -1.122189 -1.122189 0.663247 -1.123564 2.578352 -1.147887 0.453344 0.037450 1.488168
std % 0.000000 0.000000 273.973356 273.973356 332.628060 332.628060 164.360223 164.360223 340.159806 340.159806 127.655961 127.655961 638.659725 124.177131 325.572264 231.843985 710.531097 52.738871 265.561445
num 0.000000 0.000000 1.266575 1.266575 1.374615 1.374615 2.232171 2.232171 0.993610 0.993610 1.515792 1.515792 1.801816 0.603023 1.304501 0.810060 1.547864 1.131696 2.425796
diff/std % NaN NaN 182.061647 182.061647 66.790417 66.790417 123.320117 123.320117 160.751317 160.751317 74.033179 74.033179 36.809938 186.321765 197.650436 141.703951 29.288381 3.309201 61.347629
pc mean num 0.000000 0.000000 1.984662 1.984662 2.086121 2.086121 1.156253 1.156253 2.306712 2.306712 2.067538 2.067538 1.535549 2.624980 2.018566 2.563737 1.755482 1.175693 1.086203
ref num 0.000000 0.000000 1.582143 1.582143 2.352541 2.352541 1.462774 1.462774 1.536450 1.536450 1.116199 1.116199 2.476270 1.816613 0.679721 1.928507 2.621895 1.317394 0.455766
diff % NaN NaN 20.281484 20.281484 12.771050 12.771050 26.509844 26.509844 33.392197 33.392197 46.013150 46.013150 61.262895 30.795166 66.326561 24.777526 49.354677 12.052520 58.040460
num 0.000000 0.000000 0.402519 0.402519 -0.266420 -0.266420 -0.306521 -0.306521 0.770262 0.770262 0.951339 0.951339 -0.940722 0.808367 1.338845 0.635231 -0.866413 -0.141701 0.630437
std % NaN NaN 41.142539 41.142539 36.277370 36.277370 110.641730 110.641730 16.840646 16.840646 66.532772 66.532772 25.330580 11.558246 22.407086 19.911869 17.613385 76.167236 124.266573
num 0.000000 0.000000 0.816540 0.816540 0.756790 0.756790 1.279299 1.279299 0.388465 0.388465 1.375590 1.375590 0.388963 0.303402 0.452302 0.510488 0.309200 0.895493 1.349787
diff/std % NaN NaN 49.295655 49.295655 35.203903 35.203903 23.960077 23.960077 198.283347 198.283347 69.158625 69.158625 241.853501 266.434593 296.007080 124.435966 280.211202 15.823759 46.706414
In [102]:
# Plot values vs. reference cases
# NOTE - experimental code, not yet consolidated and wrapped in PEMtk

paramType = 'm'  # Currently only working for m,p cases?

# Set new DataFrame including "vary" info (missing in default case)
pDict = 'dfWideTest'
# Try using existing function with extra index set...
data._setWide(indexDims = ['Fit','Type','chisqrGroup','redchiGroup','batch', 'vary'], dataWide='dfWideTest')

# plotData = data.paramPlot(dataDict = pDict, selectors={'vary':True, 'Type':'m', 'redchiGroup':selGroup}, hue = 'chisqr', backend='hv', hvType='violin', returnFlag = True, plotScatter=True, remap = 'lmMap', hRound = 12) 
# NO REMAP CASE
plotData = data.paramPlot(dataDict = pDict, selectors={'vary':True, 'Type':paramType, 'redchiGroup':selGroup}, hue = 'chisqr', 
                          backend='hv', hvType='violin', returnFlag = True, plotScatter=True, hRound=hRound)  #, remap='lmMap')   # TODO: add renderPlot or similar option here?
p1 = data.data['plots']['paramPlot']
# p2 = dataTestSub.hvplot.scatter(x='Param',y='value', marker='o', size=200, color='green')

# Plot ref params... CURRENTLY NOT IN paraPlot(), and that also expects fit data so can't reuse directly here.

dataTest = data.data['fits']['dfRef'].copy()
# data.paramPlot(dataDict = 'dfRef')

# Set axis remap
# dataTest.replace({'Param':data.lmmu['lmMap']}, inplace=True)

# Subset
dataTestSub = data._subsetFromXS(selectors = {'Type':paramType}, data = dataTest)  
p2 = dataTestSub.hvplot.scatter(x='Param',y='value', marker='dash', size=500, color='red')

# p1+p2   # Overlays fail with "NotImplementedError: Iteration on Elements is not supported." Issue with plot types? FIXED - issues was non-plot return from paramPlot()!
p1*p2
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
*** Warning: found MultiIndex for DataFrame data.index - checkDims may have issues with Pandas MultiIndex, but will try anyway.
Out[102]:

Using the reconstructed matrix elements¶

The results tables are accessible directly, and there are also methods to reformat the best fit results for use in further calculations.

In [103]:
# self.paramsSummary contains the results above as Pandas Dataframe, usual Pandas methods can be applied.
data.paramsSummary['data'].describe()
Out[103]:
Param P_S_P_1_1_n1_1 P_S_P_1_n1_1_1 P_S_P_2_1_n1_1 P_S_P_2_n1_1_1 P_S_P_3_1_n1_1 P_S_P_3_n1_1_1 P_S_P_4_1_n1_1 P_S_P_4_n1_1_1 P_S_P_5_1_n1_1 P_S_P_5_n1_1_1 P_S_P_6_1_n1_1 P_S_P_6_n1_1_1 S_S_S_0_0_0_1 S_S_S_1_0_0_1 S_S_S_2_0_0_1 S_S_S_3_0_0_1 S_S_S_4_0_0_1 S_S_S_5_0_0_1 S_S_S_6_0_0_1
count 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000 24.000000
mean -0.386717 -0.386717 0.614867 0.614867 0.762848 0.762848 0.844967 0.844967 0.769134 0.769134 0.927061 0.927061 0.966828 1.223062 1.074484 1.149305 0.910885 0.980364 0.090045
std 1.518035 1.518035 1.169028 1.169028 1.086327 1.086327 1.288277 1.288277 1.041803 1.041803 1.238466 1.238466 1.086054 0.987656 0.987073 1.039141 1.000374 1.041297 1.482789
min -2.857393 -2.857393 -2.088247 -2.088247 -1.181846 -1.181846 -2.354861 -2.354861 -0.552797 -0.552797 0.040184 0.040184 -1.661460 -0.661391 -1.072057 -1.253769 -1.368626 -0.028552 -3.141593
25% -0.714348 -0.714348 0.152516 0.152516 0.145319 0.145319 0.184532 0.184532 0.106647 0.106647 0.090995 0.090995 0.412534 0.399965 0.377625 0.358271 0.335998 0.167667 0.032309
50% 0.064960 0.064960 0.418696 0.418696 0.458164 0.458164 0.536412 0.536412 0.393872 0.393872 0.358630 0.358630 1.328339 1.102212 1.246898 1.007032 1.315784 0.545255 0.132213
75% 0.408844 0.408844 1.142176 1.142176 1.509632 1.509632 0.846888 0.846888 1.653091 1.653091 1.113173 1.113173 1.739908 1.810177 1.851960 1.534766 1.545208 1.535203 0.220275
max 1.514140 1.514140 3.047565 3.047565 3.078588 3.078588 3.141593 3.141593 2.994350 2.994350 3.141593 3.141593 2.461075 2.992687 2.444431 3.132916 2.113221 3.141593 3.141593
In [104]:
# To set matrix elements from aggregate fit results, use `seetAggMatE` for Pandas
data.setAggMatE(simpleForm = True)
data.data['agg']['matEpd']
Out[104]:
Type m n p pc comp compC labels
Cont l m mu
P 1 -1 1 1.055450 0.255075 -2.857393 0.000000 -1.013112-0.295937j 0.255075+0.000000j 1,-1
1 -1 1.055450 0.255075 -2.857393 0.000000 -1.013112-0.295937j 0.255075+0.000000j 1,1
2 -1 1 0.754709 0.182394 -0.462299 1.984662 0.675487-0.336605j -0.073350+0.166995j 2,-1
1 -1 0.754709 0.182394 -0.462299 1.984662 0.675487-0.336605j -0.073350+0.166995j 2,1
3 -1 1 0.444572 0.107442 0.413259 2.086121 0.407146+0.178538j -0.052949+0.093488j 3,-1
1 -1 0.444572 0.107442 0.413259 2.086121 0.407146+0.178538j -0.052949+0.093488j 3,1
4 -1 1 0.697056 0.168461 1.358097 1.156253 0.147148+0.681347j 0.067851+0.154192j 4,-1
1 -1 0.697056 0.168461 1.358097 1.156253 0.147148+0.681347j 0.067851+0.154192j 4,1
5 -1 1 0.384742 0.092982 0.292101 2.306712 0.368444+0.110792j -0.062416+0.068920j 5,-1
1 -1 0.384742 0.092982 0.292101 2.306712 0.368444+0.110792j -0.062416+0.068920j 5,1
6 -1 1 0.365072 0.088229 1.187404 2.067538 0.136562+0.338568j -0.042047+0.077565j 6,-1
1 -1 0.365072 0.088229 1.187404 2.067538 0.136562+0.338568j -0.042047+0.077565j 6,1
S 0 0 0 1.650705 0.398933 0.282125 1.535549 1.585446+0.459551j 0.014058+0.398686j 0,0
1 0 0 1.434879 0.346774 0.485615 2.624980 1.268989+0.669734j -0.301519+0.171285j 1,0
2 0 0 1.513031 0.365661 0.400679 2.018566 1.393193+0.590149j -0.158315+0.329612j 2,0
3 0 0 1.356301 0.327783 0.349399 2.563737 1.274352+0.464306j -0.274563+0.179045j 3,0
4 0 0 1.345129 0.325083 0.217846 1.755482 1.313337+0.290719j -0.059698+0.319555j 4,0
5 0 0 0.483148 0.116765 2.145848 1.175693 -0.262774+0.405440j 0.044943+0.107769j 5,0
6 0 0 0.150956 0.036482 -0.913459 1.086203 0.092236-0.119500j 0.016995+0.032282j 6,0
In [105]:
# To set matrix elements from aggregate fit results, use `aggToXR` for Xarray
# data.aggToXR(refKey = 'orb5', returnType = 'ds', conformDims=True)   # use full ref dataset
data.aggToXR(refKey = 'subset', returnType = 'ds', conformDims=True)   # Subselected matE
Added dim Total
Added dim Targ
Added dim Total
Added dim Targ
Set XR dataset for self.data['agg']['matE']
In [106]:
data.data['agg']['matE']
Out[106]:
<xarray.Dataset>
Dimensions:  (Type: 2, LM: 21, Sym: 4, mu: 3, it: 1)
Coordinates:
  * Type     (Type) object 'comp' 'compC'
  * LM       (LM) MultiIndex
  - l        (LM) int64 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6
  - m        (LM) int64 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1
  * Sym      (Sym) MultiIndex
  - Cont     (Sym) object 'P' 'P' 'S' 'S'
  - Targ     (Sym) object 'S' 'U' 'S' 'U'
  - Total    (Sym) object 'P' 'U' 'S' 'U'
  * mu       (mu) int64 -1 0 1
  * it       (it) int64 1
    Eke      float64 8.1
    Ehv      float64 27.1
    SF       complex128 (2.8348434+2.8457598j)
Data variables:
    comp     (Type, mu, LM, Sym) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)
    compC    (Type, mu, LM, Sym) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)
    subset   (LM, Sym, mu, it) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)
xarray.Dataset
    • Type: 2
    • LM: 21
    • Sym: 4
    • mu: 3
    • it: 1
    • Type
      (Type)
      object
      'comp' 'compC'
      array(['comp', 'compC'], dtype=object)
    • LM
      (LM)
      MultiIndex
      (l, m)
      array([(0, -1), (0, 0), (0, 1), (1, -1), (1, 0), (1, 1), (2, -1), (2, 0),
             (2, 1), (3, -1), (3, 0), (3, 1), (4, -1), (4, 0), (4, 1), (5, -1),
             (5, 0), (5, 1), (6, -1), (6, 0), (6, 1)], dtype=object)
    • l
      (LM)
      int64
      0 0 0 1 1 1 2 2 ... 4 4 5 5 5 6 6 6
      array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6])
    • m
      (LM)
      int64
      -1 0 1 -1 0 1 -1 ... -1 0 1 -1 0 1
      array([-1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1, -1,  0,  1,
             -1,  0,  1])
    • Sym
      (Sym)
      MultiIndex
      (Cont, Targ, Total)
      array([('P', 'S', 'P'), ('P', 'U', 'U'), ('S', 'S', 'S'), ('S', 'U', 'U')],
            dtype=object)
    • Cont
      (Sym)
      object
      'P' 'P' 'S' 'S'
      array(['P', 'P', 'S', 'S'], dtype=object)
    • Targ
      (Sym)
      object
      'S' 'U' 'S' 'U'
      array(['S', 'U', 'S', 'U'], dtype=object)
    • Total
      (Sym)
      object
      'P' 'U' 'S' 'U'
      array(['P', 'U', 'S', 'U'], dtype=object)
    • mu
      (mu)
      int64
      -1 0 1
      array([-1,  0,  1])
    • it
      (it)
      int64
      1
      array([1])
    • Eke
      ()
      float64
      8.1
      units :
      eV
      array(8.1)
    • Ehv
      ()
      float64
      27.1
      array(27.1)
    • SF
      ()
      complex128
      (2.8348434+2.8457598j)
      array(2.8348434+2.8457598j)
    • comp
      (Type, mu, LM, Sym)
      complex128
      (nan+nanj) ... (nan+nanj)
      array([[[[        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj, -1.01311195-0.29593667j,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,  0.67548733-0.33660542j,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
      ...
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj]]]])
    • compC
      (Type, mu, LM, Sym)
      complex128
      (nan+nanj) ... (nan+nanj)
      array([[[[        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
      ...
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,  0.06785116+0.15419209j,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj, -0.06241586+0.06891994j,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj, -0.04204651+0.07756522j,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj],
               [        nan       +nanj,         nan       +nanj,
                        nan       +nanj,         nan       +nanj]]]])
    • subset
      (LM, Sym, mu, it)
      complex128
      (nan+nanj) ... (nan+nanj)
      dataType :
      matE
      file :
      N2O_valence.orb8_S_E0.1_4.0_28.1eV.out
      fileBase :
      /home/jovyan/N2O/N2O_valence/orb8_S
      fileList :
      N2O_valence.orb8_S_E0.1_4.0_28.1eV.out
      jobLabel :
      S/CAv
      array([[[[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]],
      
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]],
      
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]],
      
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]]],
      
      
             [[[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]],
      ...
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]]],
      
      
             [[[-0.01496218+0.01642596j],
               [        nan       +nanj],
               [        nan       +nanj]],
      
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]],
      
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]],
      
              [[        nan       +nanj],
               [        nan       +nanj],
               [        nan       +nanj]]]])

Density matrices¶

New (experimental) code for density matrix plots and comparison. See {numref}Sect. %s <sec:density-mat-basicsec:density-mat-basic> for discussion. Code adapted from the {{ PEMtk_docs }} MF reconstruction page, original analysis for Ref. {cite}hockett2023TopicalReviewExtracting, illustrating the $N_2$ case.

In [107]:
# import numpy as np

def unsignedPhase(da):
    """Convert to unsigned phases."""
    # Set mag, phase
    mag = da.pipe(np.abs)
    phase = da.pipe(np.angle)  # Returns np array only!
    
    # Set unsigned
    magUS = mag.pipe(np.abs)
#     phaseUS = phase.pipe(np.abs)  
    phaseUS = np.abs(phase)
    
    # Set complex
    compFixed = magUS * np.exp(1j* phaseUS)
#     return mag,phase
    return compFixed
In [111]:
# v2 - as v1, but differences for unsigned phase case & fix labels
# 26/07/22: messy but working. Some labelling tricks to push back into matPlot() routine

# Import routines
from epsproc.calc import density


# Compose density matrix

# Set dimensions/state vector/representation
# These must be in original data, but will be restacked as necessary to define the effective basis space.
denDims = ['LM', 'mu']
selDims = {'Type':'L'}
pTypes=['r','i']
thres = 1e-2    # 0.2 # Threshold out l>3 terms if using full 'orb5' set.
normME = False
normDen = 'max'
usPhase = True # Use unsigned phases?

# Calculate - Ref case
# matE = data.data['subset']['matE']
# Set data from master class
# k = 'orb5'  # N2 orb5 (SG) dataset
k = 'subset'
matE = data.data[k]['matE']
if normME:
    matE = matE/matE.max()

if usPhase:
    matE = unsignedPhase(matE)
    
daOut, *_ = density.densityCalc(matE, denDims = denDims, selDims = selDims, thres = thres)  # OK

if normDen=='max':
    daOut = daOut/daOut.max()
elif normDen=='trace':
    daOut = daOut/(daOut.sum('Sym').pipe(np.trace)**2)  # Need sym sum here to get 2D trace
    
# daPlot = density.matPlot(daOut.sum('Sym'))
daPlot = density.matPlot(daOut.sum('Sym'), pTypes=pTypes)

# Retrieved
matE = data.data['agg']['matE']['compC']
selDims = {'Type':'compC'}  # For stacked DS case need to set selDims again here to avoid null data selection below.
if normME:
    matE = matE/matE.max()
    
if usPhase:
    matE = unsignedPhase(matE)
    
daOut2, *_ = density.densityCalc(matE, denDims = denDims, selDims = selDims, thres = thres)  # OK

if normDen=='max':
    daOut2 = daOut2/daOut2.max()
elif normDen=='trace':
    daOut2 = daOut2/(daOut2.sum('Sym').pipe(np.trace)**2)
    
daPlot2 = density.matPlot(daOut2.sum('Sym'), pTypes=pTypes)   #.sel(Eke=slice(0.5,1.5,1)))


# Compute difference
if usPhase:
    daDiff = unsignedPhase(daOut.sum('Sym')) - unsignedPhase(daOut2.sum('Sym'))

else:
    daDiff = daOut.sum('Sym') - daOut2.sum('Sym')

daDiff.name = 'Difference'
daPlotDiff = density.matPlot(daDiff, pTypes=pTypes)



#******** Plot
daLayout = (daPlot.redim(pType='Component').layout('Component').relabel('(a) Reference density matrix (unsigned phases)') + daPlot2.opts(show_title=False).layout('pType').opts(show_title=True).relabel('(b) Reconstructed') + 
                daPlotDiff.opts(show_title=False).layout('pType').opts(show_title=True).relabel('(c) Difference')).cols(1)  
daLayout.opts(ep.plot.hvPlotters.opts.HeatMap(width=300, frame_width=300, aspect='square', tools=['hover'], colorbar=True, cmap='coolwarm'))  # .opts(show_title=False)  # .opts(title="Custom Title")  #OK



# Notes on titles... see https://holoviews.org/user_guide/Customizing_Plots.html
#
# .relabel('Test') and .opts(title="Custom Title") OK for whole row titles
#
# daPlot2.opts(show_title=False).layout('pType').opts(show_title=True).relabel('Recon')  Turns off titles per plot, then titles layout
#
# .redim() to modify individual plot group label (from dimension name) 
Set plot kdims to ['l,m,mu', 'l,m,mu_p']; pass kdims = [dim1,dim2] for more control.
Set plot kdims to ['l,m,mu', 'l,m,mu_p']; pass kdims = [dim1,dim2] for more control.
Set plot kdims to ['l,m,mu', 'l,m,mu_p']; pass kdims = [dim1,dim2] for more control.
WARNING:param.HeatMapPlot22447: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot22455: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot22484: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot22491: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot22519: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
WARNING:param.HeatMapPlot22526: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
Out[111]:

If the reconstruction is good, the differences (fidelity) should be on the order of the experimental noise level/reconstruction uncertainty, around 10% in the case studies herein.

N2O 1000 fits test case (with ~10% noise)

  • Seem to get good results for some values, particularly for imaginary terms.
  • Poor results for many real (magnitude) terms, esp. odd-L, likely because this is not defined in symmetrized dataset.
    • The main feature in the theory results around L=2,0 looks OK, but then get similar mag terms for (1,0), (3,0) and (4,0) terms, appearing as grid pattern in recon results.
    • TODO: try symmetrizing theory results to compare properly?
    • TODO: try allowing large L terms to fall off (e.g. Wigner threshold law) in fitting.

Plot MF PADs¶

Routines as per https://pemtk.readthedocs.io/en/latest/topical_review_case_study/MFPAD_replotting_from_file_190722-dist.html - currently not working. Seems to be some difference in dim stacking/assignment now...? Might be Python/Xarray version change, or PEMtk/ePSproc implementation.

In [112]:
dataIn = data.data['agg']['matE'].copy()

# # Restack for MFPAD plotter
# from epsproc.util.listFuncs import dataTypesList

# refDims = dataTypesList()
# refDims = refDims['matE']['def']
# dataStacked = dataIn.stack(refDims(sType='sDict'))
# dataStacked

# # Style 1: if full ref dataset included (with Eke dim)
# # Create empty ePSbase class instance, and set data
# # Can then use existing  padPlot() routine for all data
# from epsproc.classes.base import ePSbase
# dataTest = ePSbase(verbose = 1)

# aList = [i for i in dataIn.data_vars]  # List of arrays

# # Loop version & propagate attrs
# dataType = 'matE'
# for item in aList:
#     if item.startswith('orb'):
#         selType='L'
#     else:
#         selType = item
        
#     dataTest.data[item] = {dataType : dataIn[item].sel({'Type':selType})}

# Style 2: single Eke dim case
# Create empty ePSbase class instance, and set data
# Can then use existing  padPlot() routine for all data
from epsproc.classes.base import ePSbase
dataTest = ePSbase(verbose = 1)

aList = [i for i in dataIn.data_vars]  # List of arrays

# Loop version & propagate attrs
dataType = 'matE'
for item in aList:
    if item.startswith('sub'):
        dataTest.data[item] = {dataType : dataIn[item]}
    else:
        selType = item
        dataTest.data[item] = {dataType : dataIn[item].sel({'Type':selType})}
        
    # Push singleton Eke value to dim for plotter
    dataTest.data[item][dataType] = dataTest.data[item][dataType].expand_dims('Eke')
    
In [113]:
aList
Out[113]:
['comp', 'compC', 'subset']
In [117]:
pKey
Out[117]:
'comp'
In [123]:
# Set polarization geoms from Euler angles

# import numpy as np
# import epsproc as ep

# Set Euler angs to include diagonal pol case
pRot = [0, 0, np.pi/2, 0]
tRot = [0, np.pi/2, np.pi/2, np.pi/4]
cRot = [0, 0, 0, 0]
labels = ['z','x','y', 'd']
eulerAngs = np.array([labels, pRot, tRot, cRot]).T   # List form to use later, rows per set of angles


# Should also use MFBLM function below instead of numeric version?
# Numeric version is handy for direct surface and difference case.
R = ep.setPolGeoms(eulerAngs = eulerAngs)
# R

# Basic version - working, but get separate plots per set.
# UPDATE: use this to generate all raw figures, then restack plotly objects below

# Erange=[1,2,1]  # Set for a single E-point
Erange=[Eke-0.5,Eke+0.5,1]

# Comparison and diff
pKey = [i for i in dataIn.data_vars if i!='comp']  # List of arrays
dataTest.mfpadNumeric(keys=pKey, R = R)   # Compute MFPADs for each set of matrix elements using numerical routine
# ep.mfpad

dataTest.data['diff'] = {'TX': dataTest.data['subset']['TX'].sum('Sym')-dataTest.data['compC']['TX'].sum('Sym')}  # Add sum over sym to force matching dims
pKey.extend(['diff'])

# Plot
print(f"\n*** Plotting for keys = {pKey}, one per row ***\n")  # Note plot labels could do with some work!
dataTest.padPlot(keys=pKey, Erange=Erange, backend='pl',returnFlag=True, plotFlag=True, facetDims=['Labels',None])  # WORKING: facetDims=['Eke',None])  # DEFAULT FAILING: facetDims=['Labels', 'Eke']) # Generate plotly polar surf plots for each dataset
*** Plotting for keys = ['compC', 'subset', 'diff'], one per row ***

Summing over dims: {'Sym'}
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[compC][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.4393266707109916+0.9417986670190246j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['compC']['plots']['TX']['polar']
Summing over dims: {'Sym'}
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[subset][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.6394780210267588+0.5477229755385662j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['subset']['plots']['TX']['polar']
Summing over dims: set()
Found additional dims {'Eke'}, summing to reduce for plot. Pass selDims to avoid.
Plotting from self.data[diff][TX], facetDims=['Labels', None], pType=a2 with backend=pl.
*** WARNING: plot dataset has min value < 0, min = (-0.7282983645421228+0.4263031250551567j). This may be unphysical and/or result in plotting issues.
Set plot to self.data['diff']['plots']['TX']['polar']

Versions¶

Code¶

In [130]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [131]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews', 'pandas'])
Out[131]:
Thu Apr 11 17:15:26 2024 EDT
OS Linux CPU(s) 56 Machine x86_64 Architecture 64bit
RAM 78.6 GiB Environment Jupyter
Python 3.10.11 | packaged by conda-forge | (main, May 10 2023, 18:58:44) [GCC 11.3.0]
epsproc 1.3.2.dev0 xarray 2022.3.0 jupyter Version unknown holoviews 1.18.3
pandas 1.5.3 numpy 1.23.5 scipy 1.10.1 IPython 8.13.2
matplotlib 3.5.3 scooby 0.9.2
In [132]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
* 3d-AFPAD-dev
6fab386afe5791abfe5bf3c1c289c8c92054277f
In [133]:
# Check current remote commits
!git ls-remote --heads https://github.com/phockett/ePSproc
6fab386afe5791abfe5bf3c1c289c8c92054277f	refs/heads/3d-AFPAD-dev
226759a58ebe96bf1a6df60ccd5efb0c449af3a7	refs/heads/dependabot/pip/notes/envs/envs-versioned/pillow-10.3.0
7e4270370d66df44c334675ac487c87d702408da	refs/heads/dev
1c0b8fd409648f07c85f4f20628b5ea7627e0c4e	refs/heads/master
69cd89ce5bc0ad6d465a4bd8df6fba15d3fd1aee	refs/heads/numba-tests
ea30878c842f09d525fbf39fa269fa2302a13b57	refs/heads/revert-9-master
baf0be0c962e8ab3c3df57c8f70f0e939f99cbd7	refs/heads/testDev